# 09/10/2022 ;
# plots treatment effects from multilevel ordinal models (Figure 3) ;
library(broom)
library(dplyr)
library(ggplot2)
library(ggthemes)
library(ordinal)
library(stringr)
library(tidyverse)

theme_set(theme_tufte())

#-------------------------------------------------------------------------------
# LOAD ORDINAL LOGIT RESULTS DATASET ;
#-------------------------------------------------------------------------------
rm(list=ls())
load("TREAT_ORDINAL.RData")

# PAY ;
coef_pay_T1     <- data.frame(coef(rc_pay)$cntry[,,1])
coef_pay_T3     <- data.frame(coef(rc_pay)$cntry[,,2])
coef_pay_T4     <- data.frame(coef(rc_pay)$cntry[,,3])
coef_pay_int1   <- data.frame(coef(rc_pay)$cntry[,,4])
coef_pay_int2   <- data.frame(coef(rc_pay)$cntry[,,5])
coef_pay_int3   <- data.frame(coef(rc_pay)$cntry[,,6])
coef_pay           <- rbind (coef_pay_T1, coef_pay_T3, coef_pay_T4)
coef_pay$vignette  <- c(rep("T1: 20-25 Jahre",23), rep("T3: 50-59 Jahre",23), rep("T4: alleinerziehend",23))
coef_pay$cntry     <- substr(rownames(coef_pay),1,3)
colnames(coef_pay) <- c("ppooled","ppooled_error","ppooled_min","ppooled_max","vignette","cntry")

# EDU ;
coef_edu_T1     <- data.frame(coef(rc_edu)$cntry[,,1])
coef_edu_T3     <- data.frame(coef(rc_edu)$cntry[,,2])
coef_edu_T4     <- data.frame(coef(rc_edu)$cntry[,,3])
coef_edu_int1   <- data.frame(coef(rc_edu)$cntry[,,4])
coef_edu_int2   <- data.frame(coef(rc_edu)$cntry[,,5])
coef_edu_int3   <- data.frame(coef(rc_edu)$cntry[,,6])
coef_edu           <- rbind (coef_edu_T1, coef_edu_T3, coef_edu_T4)
coef_edu$vignette  <- c(rep("T1: 20-25 Jahre",23), rep("T3: 50-59 Jahre",23), rep("T4: alleinerziehend",23))
coef_edu$cntry     <- substr(rownames(coef_edu),1,3)
colnames(coef_edu) <- c("ppooled","ppooled_error","ppooled_min","ppooled_max","vignette","cntry")

# UNP ;
coef_unp_T1     <- data.frame(coef(rc_unp)$cntry[,,1])
coef_unp_T3     <- data.frame(coef(rc_unp)$cntry[,,2])
coef_unp_T4     <- data.frame(coef(rc_unp)$cntry[,,3])
coef_unp_int1   <- data.frame(coef(rc_unp)$cntry[,,4])
coef_unp_int2   <- data.frame(coef(rc_unp)$cntry[,,5])
coef_unp_int3   <- data.frame(coef(rc_unp)$cntry[,,6])
coef_unp        <- rbind (coef_unp_T1, coef_unp_T3, coef_unp_T4)
coef_unp$vignette  <- c(rep("T1: 20-25 Jahre",23), rep("T3: 50-59 Jahre",23), rep("T4: alleinerziehend",23))
coef_unp$cntry     <- substr(rownames(coef_unp),1,3)
colnames(coef_unp) <- c("ppooled","ppooled_error","ppooled_min","ppooled_max","vignette","cntry")

#-------------------------------------------------------------------------------
# FIGURE 3 ;
#-------------------------------------------------------------------------------

coef_pay$situation <- "S1: geringere Bezahlung"
coef_edu$situation <- "S2: geringere Qualifikation"
coef_unp$situation <- "S3: gemeinnützige Arbeit"

ppooled_coef <- rbind(coef_pay, coef_edu, coef_unp)

fig_treat <- ggplot(ppooled_coef, aes(y=vignette, x=ppooled, colour=situation)) +
  geom_point(position=position_dodge(0.5)) + 
  geom_errorbarh(aes(xmin=ppooled_min,xmax=ppooled_max, height=0), linetype="solid", size=.3, position=position_dodge(0.5)) + 
  geom_vline(xintercept=0, linetype=2) + 
  theme(legend.position="bottom", legend.title=element_blank()) +
  scale_color_grey() +
  xlab("") + 
  ylab("") +
  facet_wrap(~ cntry) + 
  scale_y_discrete(limits=c("T4: alleinerziehend", "T3: 50-59 Jahre", "T1: 20-25 Jahre") )

ggsave("FIG_3.png", plot=last_plot(), width=10, height=6)

